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The shear stress relaxation modulus G{t) may be determined from the shear stress rit) after 
switching on a tiny step strain 7 or by inverse Fourier transformation of the storage modulus G' (uj) 
or the loss modulus G" {u) obtained in a standard oscillatory shear experiment at angular frequency 
UJ. It is widely assumed that G(t) is equivalent in general to the equilibrium stress autocorrelation 
function G{t) = 13V{St (t) St(0)) which may be readily computed in computer simulations {/3 being 
the inverse temperature and V the volume). Focusing on isotropic solids formed by permanent spring 
networks we show theoretically by means of the fluctuation-dissipation theorem and computationally 
by molecular dynamics simulation that in general G{t) = Geq_ + G{t) for t > 0 with Geq_ being the 
static equilibrium shear modulus. A similar relation holds for G'{uj). G{t) and G{t) must thus 
become different for a solid body and it is impossible to obtain Geq directly from G{t). 


I. INTRODUCTION 

Background. A central rheological property charac¬ 
terizing the linear response of liquids and solid elastic 
bodies is the shear relaxation modulus G{t) [THS]. As¬ 
suming for simplicity an isotropic system, the shear re¬ 
laxation modulus G{t) = ST{t)/j may be obtained from 
the stress increment ST{t) = {f{t) — r( 0 “)) for t > 0 af¬ 
ter a small step strain with jyj <C 1 has been imposed at 
time t = 0 as sketched in panel (a) of Fig. The instan¬ 
taneous shear stress r(t) at time t may be determined ex¬ 
perimentally by probing the forces acting on the walls of 
the shear cell and in a numerical study, as shown in panel 
(b) for a sheared periodic simulation box, from the im¬ 
posed model Hamiltonian and the particle positions and 
momenta mzi. It is well known that the components 
of the Fourier transformed relaxation modulus G(t), the 
storage modulus G' {uj) and the loss modulus G"{uj), are 
directly measurable in an oscillatory shear strain exper¬ 
iment [T] as shown in panel (d) and panel (e). Using 
either G{t) or G'{uj) one obtains in the static limit (bold 
horizontal solid lines) the shear modulus [U [ 8 ] 

Geq = lim G{t) = lim G'{uj). ( 1 ) 

t—yoo UJ —^0 

The shear modulus is an important order parameter char¬ 
acterizing the transition from the liquid/sol (Geq = 0 ) to 
the solid/gel state (Geq > 0 ) where the particle permuta¬ 
tion symmetry of the liquid state is lost on the time win¬ 
dow probed mm- Examples of current interest for the 
determination of Geq include crystalline solids [TO] , glass 
forming liquids and amorphous solids mHis], colloidal 
gels m, polymeric networks UlllHHSn], hyperbranched 
polymer chains with sticky end-groups [31] or bridged 
equilibrium networks of telechelic polymers [32] . 
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FIG. 1: Sketch of notations and addressed problem: (a) 
experimental setup, (b) plain shear with periodic bound¬ 
ary conditions, (c) sketch of key equation, Eq. ([^, with 
crp = /3V{St‘^) characterizing the stress fluctuations at 7 = 0, 
(d) imposed sinusoidal shear strain 7 (t) (solid line) and mea¬ 
sured shear stress T{t) (dashed line), (e) comparison of G'{uj) 
(dash-dotted line) and storage modulus G'{uj), Eq. (j^. 


Key issue. Surprisingly, it is widely assumed [TT] [ 20 l 
|25j [28l [ 22 ] that in the linear response limit (7 ^ 0 ) 
the stress relaxation modulus G{t) must become equiv¬ 
alent to the stress autocorrelation function [4] G{t) = 
13V {ST{t)ST{0)) computed in the NV 7 T-ensemble at im¬ 
posed particle number N, volume U, shear strain 7 and 
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temperature T {f3 = \/k^T denoting the inverse temper¬ 
ature) [33]. Since G(t) = C{t) is assumed to hold, Geq 
is supposed to be measurable from some transient “finite 
frozen-in amplitude” of C{t) [20|. We call this belief “sur¬ 
prising” since, obviously, C{t) ^0 in the thermodynamic 
limit [34| for large times t 0 with 0 characterizing 
the typical relaxation time of a shear stress fluctuation 
(properly defined in Sec. VB). At variance to this belief 
we shall show by inspection of the fluctuation-dissipation 
theorem (FDT) [34| that more generally 


G{t) — Geq = C{t) for t > 0 ( 2 ) 

and G{t) = 0 for t < 0 as sketched in panel (c). Two 
immediate consequences of Eq. (§ are 

1 . that G{t) only becomes equivalent to G{t) for t > 0 
in the liquid limit where (trivially) Geq = 0 and 


II. THEORETICAL CONSIDERATIONS 
A. Reminder of some static properties 

Affine canonical displacements. Let us consider first 
an infinitessimal pure shear strain 7 of the periodic sim¬ 
ulation box in the (x,^)-plane as sketched in panel (b) 
of Fig. 12 We assume that not only the box shape is 
changed but that the particle positions r (using the prin¬ 
cipal box convention) follow the macroscopic constraint 
in an affine manner according to 

rx ^Vy for I 7 I < 1 . (5) 

Albeit not strictly necessary for the demonstration of the 
key relations we focus on, we shall assume, moreover, 
that this shear transformation is also canonical [51135]. 
This implies that the momenta must transform as 


2 . that the shear modulus Geq is only probed by G{t) 
on time scales t ^ 0 where G{t) must vanish. 

In principle, it is thus impossible to obtain the static shear 
modulus Geq of an elastic body only from G(t). We shall 
show that a similar relation 

POO 

G'(cj) — Geq = G'{uj)= / G{t) sm{ujt)d{ujt) (3) 

Jo 

pOO 

G'\uj) = G"{uj)= / G(t) cos(co’t)d(ci;t) (4) 
Jo 


Px^Px- IPy for I 7 I < 1. ( 6 ) 

All other coordinates of the positions and momenta re¬ 
main unchanged in the presented case [36] . We emphasize 
the negative sign in Eq. ^ which assures that Liouville’s 
theorem is obeyed [35] . 

Instantaneous shear stress and affine shear elasticity. 
Let liffi) denote the system Hamiltonian of a given state 
s written as a function of the shear strain 7 . The instan¬ 
taneous shear stress f and the instantaneous affine shear 
elasticity fiA may be defined as the expansion coefficients 
associated to the energy change 


holds in the angular frequency domain as sketched in 
panel (e) of Fig. [^for G'{uj). We refer below to Eqs. (§ 

as the “key relations”. 

Outline. The presented work is closely related to the 
recent paper [26] where we focus on the difference of 
static fluctuations and autocorrelation functions in con¬ 
jugated ensembles and where we also discuss briefly tran¬ 
sient self-assembled networks. The present paper pro¬ 
vides complementary informations focusing on perma¬ 
nent elastic networks in the NVyT-ensemble and on the 
response to an oscillatory shear strain 7 (t) = 7 osin(u;t). 
We begin by reminding in Sec. II A| the “affine” and 
“stress fluctuation” contributions /ta and erp to the equi¬ 
librium shear modulus Geq = Pa — and demonstrate 
then the key relations theoretically. In Sec. m we define 
our two-dimensional spring model and give some algo¬ 
rithmic details. The construction of the network and 
some properties of its athermal ground state are pre¬ 
sented in Sec. jlVI Our computational results for one finite 
temperature are given in Sec. |V| Some static properties 
are summarized in Sec. |V A We illustrate numerically in 


Sec. VB that Eq. ([^ holds. We focus in Sec. VC 


on 

the storage and loss moduli G'{uj) and G"{(jj) computed 
directly by applying a sinusoidal shear strain and com¬ 
pare this result to the Fourier transformations of G(t) 
and G(t). We summarize our work in Sec. VI where we 
briefly comment on the generalization of the key relations 
to linear response functions of other intensive variables. 


6V,ffi)/V = fy -h for ItI < 1 (7) 

with 7 = 0 being the reference, i.e. m 

T = 7i'(7)/V|^=o and ( 8 ) 

PA ^ ( 9 ) 

With Hid (7) and Hex (7) being the standard ideal kinetic 
and the (conservative) excess interaction contributions to 
the total Hamiltonian iiffi) = Hid (7) + Hex (7), this im¬ 
plies similar relations for the corresponding contributions 
Tid and fex to f = fid + t7x and for the contributions /iA,id 
and /iA,ex to Pa = Pam + AA,ex- Tho ideal contributions 
Tid and Pam are then due to the change of Hid (7) im¬ 
posed by the momentum transform Eq. (|^, the excess 
contributions fex and /iA,ex due the change of the Hex (7) 
imposed by the strained particle positions, Eq. (§. 

Ideal contributions. Using Eq. ® and p. = miv^ with 
nii being the mass of particle i ana its velocity, the 
kinetic energy of the strained system becomes Hid (7) = 
Y.i - ')ViM + 7s/ • • •]/2- This implies that 

1 ^ 

fid = -y'^miVi^xVi,y and (10) 

i=l 

1 ^ 

MA,id = 71 ) 
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for the ideal contributions to the shear stress and the 
affine shear elastiticy. Note that the minus sign for the 
shear stress is due to the minus sign in Eq. (H required 
for a canonical transformation. 

Excess contributions. We focus below on pairwise ad¬ 
ditive excess energies Eex = with u{r) being a 

pair potential and where the running index I labels the 
interaction between two particles i < j. (The poten¬ 
tial u{r) may in addition explicitly depend on the in¬ 
teraction /.) Due to Eq. a squared particle distance 
. becomes = (x + jy)‘^ + . 

Straightforward application of the chain rule [ 21 ] shows 
for the excess contributions to f and jua that 

Tex = 

I 

MA.ex = 

I 


with rii = Vi/ri being the normalized distance vector 
r = Vj — between the particles i and j. Interest¬ 
ingly, Eq. is strictly identical to the corresponding 
off-diagonal term of the Kirkwood stress tensor [4]. 

Thermodynamic averages. Let us assume an isotropic 
elastic body at imposed particle number N^ constant vol¬ 
ume E, shear strain 7 and mean temperature T (NV 7 T- 
ensemble). Note that the (intensive) shear strain 7 
corresponds thermodynamically to an extensive variable 
X = 7 E. We write F{^) for the free energy and 


Z{V = exp(-/ 3 F( 7 )) = y]exp(-/?'H( 7 )) (14) 

S 


for the corresponding partition function with being 
the sum over all accessible system states. Thermal av¬ 
erages are given by (...) = ... exp(—;d'H( 7 ))/Z( 7 ). 

Interestingly, the definition of the instantaneous shear 
stress f given above, Eq. Q, is consistent with the ther¬ 
modynamic mean shear stress r = (f), while the average 
affine shear elasticity //a = (Aa) only corresponds to an 
upper bound to the shear modulus Geq. To see this let us 
first note for convenience that 


dlogjZj'f)) 

dj 

log(Z( 7 )) 

972 

and 


Z'il) 

Z{i) 

Z"h) fz'ivy 

^(7) V ^(7) J 


(15) 

(16) 


Z'h) = -^/3?i'(7)e-^«(^) (17) 

s 

Z"iV = ( 18 ) 


It then follows using Eq. (15) and Eq. that the mean 
shear stress is indeed m 


dF{X) 


dx 


x=o 


(19) 


= = if). ( 20 ) 


Using Eq. (16) and Eq. ( [T^ one verifies also that 
d'^F{X) 


Ge. = U 


( 21 ) 


x=o 


= Ma — ctf with cTp = PV {Sr‘^) ( 22 ) 

characterizing the shear-stress fluctuations at 7 = 0 . 
Since crp > 0, /iA is only an upper bound to Geq. 
We emphasize that Eq. ( [ 2 ^ is a special case of the 
general stress-fluctuation relations for elastic moduli 
[39 1 140]. It provides a computational convenient 
method to obtain Geq for systems in the NVyT-ensemble 
used in many recent numerical studies [EHISEIIESIIIS]. 
We have demonstrated Eq. ( 22 ) without introducing a lo¬ 
cal displacement field as in Ref. [39] . We note en passant 
that the averages r and fiA are “simple averages”, i.e. 
no fluctuations, and are thus identical for any ensemble 
given that the same state point is sampled BED- 

Lebowitz-Percus- Verlet transform. Equation ( 22 ) can 


alternatively be obtained from the general transforma¬ 
tion relation for a fluctuation (SASB) of two observables 
A and B due to Lebowitz, Percus and Verlet im 


sAss'j 


= (6A6B 


dm d{A) d{B) 
X dX dipl)d{pl) 


(23) 


with X = denoting again the extensive variable and 
I = r the conjugated intensive variable [42]. This gives 


(JfI^ = CTfI^ + Geq (24) 

i.e. the thermodynamic shear modulus Geq compares the 
shear stress fluctuations in the conjugated ensembles at 
constant mean shear stress r and imposed shear strain 7 . 
The latter formula can be made more useful for compu¬ 
tational studies by rewriting the shear stress fluctuations 
(JfI^ at constant shear stress r. Note that the normal¬ 
ized weight of a state s in the NVrT-ensemble is given by 
p{j) ^ exp[—/ 3 (H( 7 ) — Eyr)]. Using the instantaneous 
shear stress f defined in Eq. Q we thus have 

P'il) = -l3V[T{-y) - t)]p(7). (25) 


Using integration by parts it is then readily seen m that 
this leads to = {'H"{'y)/V)\r = {P"{^)/V)L = /ta 
in agreement with Eq. This confirms Eq. ( [2^ [33] . 

Simplifications. Thermal averaging of Eq. (jlQj) and 
Eq. (II) implies rid = /iA,id = ^id with Pid being the ideal 
normal pressure. We note further that crp = ap^d + crF,ex 
may again be rewritten as the sum of an ideal con¬ 
tribution cTpqd = /3V {6r^fi) = Pid and an excess term 
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c^F,ex = /3V{6t^^). All ideal contributions to Geq thus 
cancel and one may rewrite Eq. (22) as 


Geq /^A,ex <^F,( 


(26) 


used for the determination of the correlation functions 
[4l l4n[4^ . Generalizing the relation between the static 
fluctuations Eq. (24) into the time domain, this yields 
immediately 


Since an ideal gas must have a vanishing shear modulus, 
this simplification is, of course, expected and could have 
been used from the start [21] . Using Eq. (10) and the fact 
that in an isotropic system all coordinates are equivalent 
it is seen that the average shear elasticity reduces to m 


MA,ex = MB - -Pex with 

Mb = 4^((rfu"(r0(27) 

being the Born-Lame term [40l [43] and Pex the excess 
contribution to the normal pressure P = Pid + Pex- 


B. Demonstration of Eq. @ 

Static limits. As shown by the dash-dotted line in 
panel (c) of Eig. by definition G(t) ^ crp for t ^ 0 
and G(t) ^ 0 for t ^ oo [34]. Equation ([^ thus implies 
that the relaxation modulus becomes 


G(t) — y (Jp -|- (/Ta — <^f) — /^A ^ ^ 


which is consistent with an affine canonical shear strain 
imposed at t = 0, Eq. ([^, and G(t) ^ Geq for t ^ oo as 
expected from Eq. §• 

Fluctuation Dissipation Theorem. We show next that 
Eq. ([^ must hold for all times. Using Boltzmann’s su¬ 
perposition principle for an arbitrary strain history 7 (t) 
[1] the shear stress becomes [3] 

= [ dg G{t — s) (29) 

J-oc ds 

= G{t - s)j{s)\*_^ - J 


using integration by parts in the second step. Introduc¬ 
ing the “after-effect function” or “dynamic response func¬ 
tion” [37] x(t) = —G'{t) = G'{—t) this may be rewritten 
for a step strain imposed at t = 0 as 


/ oo 

x{s) ds 


(30) 


being consistent with the expected G{t) Geq for 
t ^ oo. Since according to the EDT as formulated 
by Eq. (7.6.13) of Ref. [34], the after-effect function is 
x{t) = —G'{t)^ this yields the claimed relation Eq. (©• 
Alternative demonstration. It is of some importance 
that Eq. (§ may be also obtained from the general trans¬ 
formation relation Eq. (23) with A = r(t) and P = r(0). 
It is assumed here that the shear-barostat imposing a 
mean shear stress r = 0 is sufficiently slow such that 
the system trajectory is not altered over the time scales 


+ ^eq- (31) 


As before for the static shear stress fluctuations crpl^ one 
can show that 



dx 


T 


G{t) for t > 0 


(32) 


where we have reexpressed in the first step [f(t; 7 ) — 
r][r( 0 ; 7 ) — t]p{x) using Eq. ( [^ . In the second step 
we have used that within linear response G(t) does not 
depend on 7 . Using Eq. (|32|) and G{t) = G(t) |^, Eq. (31) 


implies again Eq. SET 



C. Oscillatory shear 


Experimentally, the relaxation modulus G(t) is, of 
course, commonly sampled in a linear viscoelastic mea¬ 
surement mm using an oscillatory shear imposing, e.g., a 
sinusoidal shear strain x{t) = 7 o sin(co’t) of amplitude 70 
and angular frequency uj as shown in panel (d) of Eig. 
This implies an average shear stress 

r{t) = Tref + 7 o (G'(cj) sin(u;t) + G"{uj) cos{ujt)) (33) 


with Tref being the (not necessarily vanishing) reference 
shear stress at 7 = 0. The storage modulus G'{uj) and 
the loss modulus G"{uj) may thus be determined as the 
Eourier coefficients 

G'iuj) = [ sin(u;t)^-^dt (34) 

pTuj Jo To 

G"{uj) = ^ / cos(cjt)^dt (35) 

pTco Jo 70 

with = l/27ruj being the period of the oscillation. 
Averages over instantaneous shear stresses T(t) are per¬ 
formed over all periods p sampled. As stated, e.g., by 
Eq. (7.149) and Eq. (7.150) of Ref. [1] both moduli are 
on the other side quite generally given by the Eourier-Sine 
and Eourier-Cosine transforms of G(t) 


G'{lj) — Geq 
G"{uj) 



Geq) sin(u;t)dt, (36) 
Geq) cos{ujt)dt. (37) 


The latter two relations together with Eq. ([^ imply the 
key relations Eq. and Eq. (|^ announced in the Intro¬ 
duction. We shall pay special attention to the low-cj and 


high-cj limits of G'{uj) and G"{uj) at the end of Sec. V C 
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FIG. 2: Snapshot of a small subvolume of linear length 10 
containing about 100 verticies of the elastic network consid¬ 
ered in this work assuming Eq. (38). The total periodic box of 
linear length L ^ 102.3 contains N — 10^ verticies and 9956 
harmonic springs. The light grey circles indicate the positions 
and diameters of the quenched polydisperse LJ bead system 
which was used for the construction of the elastic network as 
described in Sec. El The lines represent the quenched forces 
of the athermal (T = 0) reference configuration. The dark 
(black) lines indicate repulsive forces between the verticies, 
while the light (red) lines represent tensile forces. The line 
width is proportional to the tension/repulsion. 


III. MODEL AND ALGORITHMIC DETAILS 

Hamiltonian. To illustrate our key relations we 
present in Sec. |V| numerical data obtained using a pe¬ 
riodic two-dimensional (d = 2 ) network of harmonic 
springs. The model Hamiltonian is given by the sum 
^ = ^id + ^ex of a kinetic energy contribution Hid = 
(assuming a monodisperse mass m) and an 
excess potential 

Tiex = ^Ui{ri) with ui{r) = ]^Ki (r - Rif (38) 

i 

where Ki denotes the spring constant, Ri the reference 
length and ri = |r^ — r^ | the length of spring 1. The 
sum runs over all springs I connecting pairs of beads i 
and j at positions and r^. A small subvolume of the 
network considered is represented in Fig. An experi¬ 
mentally relevant example for such a permanent network 
is provided by endlinked or vulcanized polymer networks 
[U |28l [29] . Since the network topology is by construction 
permanently fixed, the shear response G{t) must become 
finite for t ^ oo for all temperatures at variance to sys¬ 
tems with plastic rearrangements as considered, e.g., in 
Ref. [To]. Note that the particle mass m and Boltzmann’s 
constant /cp are set to unity. Lennard-Jones (LJ) units 
[4] are assumed throughout this paper. 


Computational methods, parameters and observables. 
The construction and the characterization of the refer¬ 


ence network at zero temperature is presented in Sec. IV 


As discussed in Sec. |V| this network is then investi¬ 
gated numerically by means of standard molecular dy¬ 
namics (MD) simulation [HIS] at constant particle num¬ 
ber N = 10^, box volume V ~ 102.3^ and a small, but 
finite mean temperature T = 0.001. Newton’s equations 
are integrated using a velocity-Verlet algorithm with a 
tiny time step dtuB = 10“^. The temperature T is fixed 
using a Langevin thermostat with a friction constant 
ranging from C = 0.01 up to C = 5 as specified below. An 
important property sampled in this study is the instan¬ 
taneous stress tensor HE?] 

m ^ 1 

<7a/3 = -y'^Vi^aVi,^ + (39) 


i=l 


V 


where the Greek indices a and fi denote the d = 2 spa¬ 
tial directions x and y, Vi^a the ^-component of the ve¬ 
locity of particle i and ni^a the ^-component of the nor¬ 
malized vector r^. The first term in Eq. (39) gives the 
ideal gas contribution, the second term corresponds to 
the Kirkwood excess stress for pair interactions. The 
trace of the stress tensor yields the instantaneous nor¬ 
mal pressure P = —{cFxx + ^yy)/^^ the off-diagonal el¬ 
ements CTxy = ^yx ^rc consistcnt with the shear stress 


^ = 'fid + fex, Eqs. ( 1 Q| 12 ), obtained in Sec. EH The 
shear strain 7 is set to zero for the computation of the 
autocorrelation function C{t) = fiV ((T(t)^) — A 

tiny step strain 7 <C 1 is imposed at t = 0 in order to 
compute G(t), as discussed in Sec. VB[ and a sinusoidal 
strain ^(t) = 70 sm{ujt) to measure directly the storage 
and loss moduli G'{uj) and G"{uj) from the sine and co¬ 
sine transfo rms of the instantaneous shear stress f(t), 
Eqs. ( |34|35| ). As already described in the first paragraph 
of Sec. |II A[ we perform in both cases affine, canonical [35] 
and (essentially) infinitessimal strain transformations of 
the b ox s hape and the particle positions and momenta, 
Eqs. ( |5|6| ). Note that the transformation of the momenta 
is, however, not important for the simulations presented 
here where we focus on low temperatures and mainly on 
high values of the Langevin friction constant. Moreover, 
as shown in Ref. [26], similar results, especially concern¬ 
ing Eq. are also obtained using Brownian dynamics 
or Monte Carlo simulations [ilS]. 


IV. REFERENCE CONFIGURATION 

Construetion of network. As explained in detail in 
Ref. [ 21 ] , our network has been constructed using the dy¬ 
namical matrix M of a polydisperse LJ bead glass com¬ 
prising V = 10^ particles. Prior to forming the network 
a LJ bead system has been quenched to T = 0 using a 
constant quenching rate and imposing a normal pressure 
P = 2. The original LJ beads are represented in Fig.j^by 
grey polydisperse circles, the permanent spring network 
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FIG. 3: Energy difference per volume e after an affine strain 
7 is applied (circles) and after the network has relaxed to the 
new ground state (squares). Main panel: Double-logarithmic 
representation of e vs. 7 . The thin solid line indicates 
Eq. (40), the bold solid line Eq. (41). Inset: Half-logarithmic 
representation of reduced energy increment (e — t 7 )/( 7 ^/ 2 ) 
vs. 7 where the residual shear stress r = 0.011614 has been 
taken off. The shear modulus Geq ~ 16.2 is obtained from 
the rescaled groundstate energy (squares). 


created from the quenched bead system by lines between 
verticies. The network has a number density p ~ 0.96 
corresponding to a linear periodic box length L ^ 102.3. 
Using Eq. (39) one determines a normal pressure P ^ 2 
and a small, but finite shear stress r ~ 0.011614 (deter¬ 
mined to high precision for reasons given below). Using 
Eq. (27) one obtains the affine shear elasticity ~ 34.3. 


By construction the total force acting on each vertex of 
the reference network vanishes albeit the repulsive and 
tensile forces transmitted along each spring do in general 
not as shown in Fig. Due to the periodic boundary 
conditions the finite residual (normal and shear) stress 
does not relax. Albeit being small, it is important to 
properly account for the finite, quenched shear stress r of 
the reference for all the correlation functions considered. 
We also note that the force network. Fig. is strongly 
inhomogeneous with zones of weak attractive links em¬ 
bedded within a strong repulsive skeleton as discussed in 
Refs. pMl[T4] . 

Affine displacements. As shown in Fig. it is possi¬ 
ble to obtain numerically the shear stress r and the affine 
shear elasticity from the energy increment per volume 
e = (Pex( 7 ) — ^ex( 0 ))/U caused by an affine shear strain 
7 (spheres). In agreement with Eq. 0 , the energy in¬ 


creases as 


e = T 7 + /ua7^/2 


(40) 


as shown by the thin solid line in the main panel of Fig.|^ 
The asymptotic limits for small and large 7 are indicated 
by, respectively, the dashed line and the dash-dotted line. 
We have used for the coefficients r and pA in Eq. (40) 



EIG. 4: Snapshot of non-affine displacement field after an 
affine shear strain with 7 = 0.01 is applied and the system 
is allowed to relax. The arrow length is proportionel to the 
non-affine particle displacements. The displacement field is 
correlated over distances much larger than the typical particle 
distance. All displacements are strictly linear in 7 . 


the values obtained using Eq. (12) and Eq. (13). This is 


merely a self-consistency check since both properties are 
actually defined assuming a virtu al affine strain transfor¬ 
mation as reminded in Sec. Infill E]. As shown in the 
inset, an accurate verification of pA over essentially the 
full 7 -range is obtained by plotting in half-logarithmic 
coordinates the reduced energy (e — /2). Eor the 

smallest 7 an even more precise value of the substracted 
residual shear stress r is required to collapse even these 
data points on the dash-dotted horizontal line. 

Non-affine displacements. The forces /. acting on the 
verticies i of an affinely strained network do not vanish 
in general. Eollowing Ref. m the network is first re¬ 
laxed by steepest descend, i.e. imposing displacements 
proportional to the force, and then by means of the con¬ 
jugate gradient method [ 6 ]. These non-affine displace¬ 
ments lower the final energy of the strained system as 
indicated by the squares in Eig.[^ As shown by the bold 
solid line in the main panel, these final energies scale as 


— T7 + Geq7^/2. 


(41) 


This scaling is similar to the affine strain energy, Eq. ( [4Q| ) , 
having the same linear term but with pA being replaced 
by the shear modulus Geq. As shown in the inset, these 
energies can thus be used to determine Geq ~ 16.2 by 
taking again into account the quenched shear stress at 
7 = 0 . It follows from Eq. (40) and Eq. (ETl) that 


due to the non-affine displacements a substantial fraction 
1 ~ Geq/P a ~ 1/2 of the affine strain energy is relaxed 
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FIG. 5: Determination of shear modulus Geq. Inset: Low¬ 


est eigenfrequencies ujp obtained by diagonalization of the 
dynamical matrix. The bold horizontal lines correspond to 
Gealp ~ 16 according to Eq. (43) fitting the eigenvalues for 
p = 3,4, 5, 6. Main panel: Fourier transformed transverse 
displacement field for a temperature T = 0.001 as a func¬ 
tion of wavevector q. The horizontal line corresponds again 
to Geq_ ~ 16. The open symbols indicate data obtained us¬ 
ing the ground state positions as reference, the filled symbols 
assume the average vertex position as reference. 



FIG. 6: Eigenvector for p = 4 (n = l,m = 0) correspond¬ 
ing to a transverse planewave with wavector q pointing in 
the horizontal direction, i.e. with displacements essentially in 
the vertical direction, and wavelength \ — L. Interestingly, 
heterogeneous deviations from the planewave solution of con¬ 
tinuum mechanics are even visible for this low eigenmode. 


for large 7 in agreement with Ref. M- A snapshot of 
these non-affine displacements is given in Fig. Note 
that the non-affine displacements are correlated over dis¬ 
tances much larger than the typical particle distance and 
one thus expects deviations from standard continuum 
mechanics if similar length scales are probed. 

Eigenstates of the referenee. Following Refs. [13 El 
the shear modulus Geq of the network at T = 0 may 
alternatively be computed from the lowest non-trivial 
eigenfrequencies ujp. (The running index p increases 
with frequency.) This is shown in the inset of Fig. 
The eigenfrequencies have been determined by di¬ 
agonalization of the dynamical matrix M by means of 
Lanczos’ method [ 6 ]- It follows from continuum elastic¬ 
ity El that the eigenfunctions must be planewaves with 
wavevectors q quantified by the boundary conditions, i.e. 
{Qx^ Qy) = (27r/I/) (n, m) with n, m = 0, ±1 ,... being two 
quantum numbers. An example for such a planewave is 
given in Fig.j^for p = 4 and the pair of quantum numbers 
(n,m) = (1,0). For the wavelength this implies 

\p = 27 r/|g| = L/^ w? (42) 

and for the eigenfrequencies of the transverse modes 

cjp = with c_l = (43) 


being the transverse wavevelocity. As shown by the hori¬ 
zontal lines, we obtain by fitting Eq. (43) to the frequen¬ 
cies for p = 3,4, 5, 6 (corresponding torr + m? = 1) that 


Geq ~ 16 and ~ 4. Interestingly, the degeneracy of 
the eigenvalues expected from continuity elasticity is al¬ 
ready lifted for p = 7, 8 , 9,10 (n^ + = 2), i.e. the box 

size L does not allow a precise determination of Geq using 
these eigenvalues. Deviations from the planewave solu¬ 
tion are even visible from the eigenvector for p = 4 pre¬ 
sented in Fig. Continuum elasticity must break down 
in any case if wavelengths of order of about the particle 
distance are probed, A ^ 1. This implies that Eq. (43) 
can at best hold up to a frequency uj = 27rc^/A ^ 25. 
We come back to this issue in Sec. IV Al and Sec. IV Cl 


V. FINITE-TEMPERATURE 
COMPUTATIONAL RESULTS 

A. Some static properties 

Introduetion. Since the temperature T = 0.001 is 
rather small, one expects all static properties such as the 
pressure P or the elastic modulus Geq to be similar to 
their ground state values. As we have checked comparing 
various methods one confirms indeed that P ^ Pex ~ 2 , 
r « Tex ~ 0.011, /iA ~ /iA,ex ~ 34, CTp ~ CrF,ex ~ 18 
and Geq ~ 16. The same applies in fact to all small 
temperatures T ^ 1. 

Displaeement eorrelations. A finite-T method for 
computing Geq, which is also of experimental relevance 
[ 20 ] . is presented in the main panel of Fig. Using 
an ensemble of 10 ^ configurations sampled over a total 
time interval 10^ with C = 5 we obtain first the dis- 
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FIG. 7: Determination of Geq ~ 16 using Eq. ( [^ as a func¬ 
tion of the sampling time t at T = 0.001. While the simple 
average /xa ~ 34 converges immediately (small filled squares), 
the stress fluctuations crp (t) are seen to increase monotonously 
to the large-t limit crp ~ 18 (bold solid line). This is shown 
for four different friction constants = 5,1, 0.1 and 0.01. The 
thin solid line compares the data for ( — 5 with the integrated 
correlation function C{t) using Eq. ( [4^ . 


placements for each vertex particle i of the network 
either by taking the position of the ground state net¬ 
work as reference for defining the displacement (open 
squares) or alternatively the average particle position 
in the ensemble (filled squares). The displacement field 
m(e) = '^i U.i^{L — Li) is then Fourier transformed accord¬ 
ing to u{q) = Yl!i=i exp(ig' • r) yJ(r)/\/N [20]. Note that 
the wavevector q must be commensurate to the square 
simulation box of linear length L, Eq. ( [4^ . The compo¬ 
nent of u{q) perpendicular to q corresponds to the trans¬ 
verse component u^{q) of the Fourier transformed dis¬ 
placement field. Using that according to the equiparti- 
tion theorem every independent elastic mode corresponds 
to an average kinetic or potential energy /cbT/ 2 , contin¬ 
uum mechanics implies that 0120 ] 

2 /= (|m±( 2 )|^)-> fcBTp/Geq forg->0 (44) 

where the average (...) is taken over the ensemble of 
stored configurations. As can be seen from Fig. 0 V 
becomes rapidly constant below q ^ 1 confirming Geq ~ 
16 as indicated by the bold horizontal line. Incidentally, 
both definitions of the displacement field yield identical 
results. Since only the second definition using the average 
particle position can normally be used in experimental 
studies, this is rather reassuring. Note also that the first 
peak for large wavevectors corresponds to a wavelength 
A = 27r/g' ~ 1 of order of the typical particle distance 
where continuum mechanics should break down. 

Stress fluctuations. While it is thus possible to get 
Geq from the displacement correlations, it is from the 
computational point of view more convenient to deter¬ 
mine the modulus using the stress-fluctuation formula. 


FIG. 8 : Stress relaxation modulus G{t) and autocorrelation 
function G(t) sampled by means of a Langevin thermostat of 
friction constant = 5. An affine step strain 7 = 0.01 is 
applied at t = 0 to determine G{t). 


Eq. (22). This is shown in Fig. where the affine 
shear elasticity jua (small filled squares) and the stress- 
fluctuation term a^it) are presented as functions of the 
sampling time t. (The notation ap without time argu¬ 
ment refers to the static thermodynamic large-t limit, 
while crF(t) indicates that this property has been deter¬ 
mined using a finite time window t.) While the sim¬ 
ple average /a a is obtained immediately, crF(t) is seen to 
increase monotonously from zero to the large-t plateau 
cTp ~ 18 (horizontal bold solid line). Friction effects are 
important only in the intermediate sampling time win¬ 
dow between t ~ 0.1 and 10 where the stress fluctuations 
strongly increase. The data obtained for Langevin fric¬ 
tion constants C = 0.1 and 0.01 are essentially identical. 
The time dependence of crF(t) for different friction con¬ 
stants is further analysed at the end of Sec. lyi 


B. Time domain of key relations 

High friction limit. We turn now to the numerical 
verification of the key relations stated in the Introduc¬ 
tion focusing first on results obtained in the time domain 
which have already been reported in Ref. [26]. The data 
presented in Fig. have been computed using a Langevin 
thermostat with one high friction constant C = This 
is done merely for presentational reasons since in this 
limit the kinetic degrees of freedom can be disregarded 
and oscillations are suppressed. The stress relaxation 
modulus G(t) given has been computed as indicated in 
panel (b) of Fig. by applying an affine canonical shear 
strain 7 = 0.01 and by averaging over 1000 runs start¬ 
ing from independent configurations at t = 0“. Due to 
the strong damping G(t) decreases monotonously from 
G( 0 +) = /iA to Geq while C{t) decays from G( 0 ) = crp 
to zero. Confirming Eq. ( 0 . only after vertically shifting 
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FIG. 9: Stress relaxation modulus G(t) obtained using a step 
strain (open symbols) and the rescaled autocorrelation func¬ 
tion Geq_ G{t) at 7 = 0 (filled symbols) for several friction 
constants Inset: Shear stress relaxation time 0{(^) measured 
using Oi/e, 0^% and 0i% defined in the main text. 


C{t) C (t) H- Geq one obtains a collapse on the directly 
computed modulus G{t). This is the most important nu¬ 
merical result of the present work. 

Friction effects and stress relaxation time. A simi¬ 
lar scaling collapse of G{t) and G{t) + Ggq has been 
also obtained for different friction constants ( as may 
be seen from Fig. As one expects, the MD data de¬ 
cays more rapidly with decreasing ( and reveals anti¬ 
correlations for the lowest ( probed. The effect of the 
friction constant may be compactly described using the 
stress relaxation time 0 shown in the inset of Fig. We 
compare here three different definitions. The definition 
G{t = 6 >i/e) = cTp/e, indicated by spheres, describes 
the time where the correlation functions first becomes 
small. This simple definition underestimates correlations 
at larger times. Alternatively, one may attempt to define 
0 by the ratio Ci/cq of the first to the zeroth moment of 
G{t) with Cn = t'^G{t)dt (not shown). Unfortunately, 
Cl is found to converge badly and we have not been able 
to obtain reliable values over the full C-range. A nu¬ 
merically well-behaved alternative is based on an exact 
relation between the autocorrelation function G{t) and 
the monotoneously increasing fluctuation term crF{t) de¬ 
termined as a function of the sampling time t. Assuming 
time-translational invariance it can be shown [26] that 


y{t) = 1 - = jf ds (1 - s/t) 

(T'f t Jo 


C{s) 

CTp 


(45) 


i.e. the dimensionless ratio y{t) characterizes the dif¬ 
ference of the zeroth and the first moment of the au¬ 
tocorrelation function integrated up to t. That this re¬ 
lation holds can be seen from the thin line presented in 
Fig. 0 We note en passant that y(f) must ultimately van¬ 
ish slowly as 1 /t since the integral over G{s) in Eq. (45) 


becomes constant [26]. (This 1/t-decay is consistent with 
the general finite-sampling time corrections for fluctua¬ 
tions [7]-) Since (JY{t) can be determined directly, one 
may define ^(0 using a fixed ratio y{0). We have pre¬ 
sented in Fig. m the constants y{0^%) = 5% (squares) 
and y{0i%) = 1% (triangles). Note that 0^% corresponds 
to the time where cr-pit) begins to saturate, while 0i% 
indicates the time where the stress-correlations become 
neglible and a^it) thus allows a good estimation of Geq 
using the stress-fluctuation formula. 


C. Oscillatory shear 


Direct determi natio n ofG'{uj) and G"{uj). As already 
stressed in Sec. IIC, the relaxation modulus G{t) is 
commonly determined experimentally by inverse Fourier 
transformation of the storage modulus G'{uj) and/or the 
loss modulus G"{uj) obtained in a linear viscoelastic mea¬ 
surement imposing an oscillatory shear mm- Motivated 
by this we perform non-equilibrium MD simulations in 
the linear response limit by imposing a sinusoidal shear 
strain of frequency uj with an amplitude 70 = 0.001. (By 
varying 70 it has been checked that all reported values 
are in the linear regime.) This is done by performing 
every time step dtuB = 10 “^ an affine canonical strain. 
Sec. |II A| We use production runs of total length pT^ 
with = 27t/u and at least p = 100 oscillation peri¬ 
ods, i.e. the computational load increases as p/uj with 
decreasing frequency. This sets a lower limit 00 « 10“^ 
for the angular frequencies we have been able to sam¬ 
ple. Two production runs are compared to rule out tran¬ 
sient behavior. The instantaneous shear stress f(t) is 
sampled each StuB- Using Eq. (34) and Eq. (j^ one 


obtains G'{uj) and G"{uj) as indicated by open triangles 
in Eig. where we focus for simplicity on one value 
C = 5 of the friction constant, and large open symbols in 
Fig.[TTJ where data for different (f are presented. 


Test of key relations. These values of G'{uj) and 
G"(uj) are compared in Eig.j^with the Sine-Eourier and 
Cosine-Eourier transforms of the shear modulus G{t) ob¬ 
tained from the step strain experiment (circles) and of 
the shear stress autocorrelation function G(t). To reduce 
the well-known (but considering our tiny time step 
surprisingly severe) numerical problems at large frequen¬ 
cies {uj ^ 10), Eilon’s quadrature method [4] is used for 
the Eourier integration. These Eourier transforms col¬ 
lapse nicely on the respective directly computed moduli 
G'{uj) and G"{uj). The crucial point is here that it is 
Geq + G'{uj) which collapses on G'{uj)^ while G'{uj) does 
not being much too small considering that Geq ~ 16, as 
seen in panel (a). We have verified that a similar data col¬ 
lapse for Geq+G'(c(;) onto G' {uj) and G"{uj) onto G" {uj) is 
also obtained for other friction constants as may be seen 
from Eig. [^ Having settled this fundamental scaling 
issue let us turn finally to the description of the general 
shape and asymptotic behavior of both moduli. 
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FIG. 10: Oscillatory shear moduli (a) G'{uj) using a log-linear 
representation and (b) G"{(ju) using a double-logarithmic rep¬ 
resentation for a large friction constant ^ = 5. We compare 
the Fourier coefficients (triangles) obtained by applying a si¬ 
nusoidal shear strain with amplitude 70 = 0.001 to the Fourier 
transformed step-strain relaxation modulus G{t) (circles) and 
shear autocorrelation function G{t) (squares). The bold solid 
lines indicate the expected asymptotic limits for small u, the 
dashed lines the large-ct; asymptotics, the thin vertical lines 
the frequencies uj ^ 25 and oj ~ 15 of the maxima of both 
moduli. The latter values coincide with the breakdown of 
continuum elasticity at a planewave frequency u = 27rc^/A 
with X ^ 1 probing the typical particle distance. Inset of (b): 
Phase angle S with maximum S ~ 20° at Ct; ~ 7 (vertical line). 


Low-frequency limit. The low-cj limit of G'{uj) and 
G"{uj) are readily obtained by expanding sin(x) in Eq. @ 
and cos(x) in Eq. As one expects [T] , one obtains to 
leading order that G'{uj) Ggq, Eq. (1), and G"{uj) 
ujcq as indicated by bold solid lines in Eig.[TQ]and Eig.pTj 
We have fitted G"{uj) for (" = 5 using the independently 
determined value cq ~ 16. 

Intermediate peaks. In the intermediate frequency 
regime one observes a striking peak at cj ~ 25 for G'{uj), 
i.e. at the planewave frequency corresponding to the typ- 


FIG. 11: Effect of Langevin friction constant ( with open sym¬ 
bols referring to data obtained using oscillatory shear with 
amplitude 70 = 0.001: (a) Storage modulus G'{uj). The fric¬ 
tion constant ( only matters in an intermediate frequency 
range between a; 0.1 and ca 10 below the peak at cj 25. 
(b) Loss modulus G"{uj). The dot-dashed line is a phe¬ 
nomenological fit for large uj and small C assuming Eq. ( [47| . 


ical particle distance, Eq. ( [4^ , and at a slightly smaller 
value UJ ^ 15 for G"{uj). The phase angle 6 characteriz¬ 
ing the ratio G" {uj)/ G' {uj) = tan(J) is presented in the 
inset (in grad). It reveals a maximum at u; ~ 7. 

High-frequency limit. In the high-frequency limit it 
becomes increasing difficult to relax the imposed affine 
strain by non-affine displacements and, hence, to dissi¬ 
pate the work done on the system. This implies that 

G' {uj) = Geq + G'{uj) (/iA ~ CTp) + (Jp = /^A (46) 

for cj ^ oc (bold dashed horizontal lines), while G"{uj) 
decays with a non-universal (see below) sharp cutoff 
(dashed and dash-dotted lines). This asymptotic behav¬ 
ior can be also understood on mathematical grounds by 
expanding Eq. © and Eq. 0 by integration by parts. 
This leads to two series-expansions in terms of even and 
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odd derivatives of C{t) taken at t = 0 for, respectively, 
G'{uj) and We remind that G{t) is an even func¬ 

tion for classical systems [4l|34] and that in a MD simula¬ 
tion G{t) must become analytic at sufficiently short times 
t 1/C) where thermostat effects are negligi¬ 
ble. (Strict non-analyticity of G{t) at t = 0 is formally 
found, e.g., for the Maxwell model or the polymer Rouse 
model [3].) Analyticity at t = 0 implies that the (con¬ 
verging) expansion of G'{uj) is dominated by its first term 
for cj ^ oo. This demonstrates Eq. (46). The effective 


odd derivatives of G{t) must, however, become negligible 
with decreasing friction C, increasing angular frequency 
uj and decreasing time step dtuB- Hence, G"{uj) must 
rapidly vanish. This point is further investigated in the 
final paragraph of this section. 

Friction effects. The friction effects presented in 
Fig. 11 are rather small for G'{uj) due to the large static 
constants Geq and fiA dominating the storage modulus. 
As one would expect, they are more marked for the loss 
modulus G''{uj), especially at low frequencies where the 
thermostat has more time to dissipate the applied me¬ 
chanical energy, but also at the large-u; cutoff. For C = 5 
one might be tempted to fit an 1/cj^-decay (dashed lines) 
which indicates an effective non-analytical behavior of 
G{t) at short times. With decreasing C the loss modulus 
is, however, better described by the exponential decay 


G"{uj) = cTp /4 X exp(—x^/2) with x = uoO (47) 


as indicated by the dash-dotted line in panel (b) of 
Fig. [TT] This phenomenological fit is the Cosine- 
Fourier transform assuming a Gaussian G{t) = 
cTp exp(—(t/^)^/2) with the parameter 0 ^ 0.1 deter¬ 
mined from the short-time (high-frequency) limit. A 
similar phenomenological fit for the large-cj behavior 
of is also obtained by assuming a Lorentzian G{t) = 
cTp / (1 + (t/^) ^). We note finally that if smaller time steps 
and larger frequencies could be computed, one expects 
even for higher friction constants a similar non-algebraic 
cutoff, albeit with a different parameter 0{(). 


VI. CONCLUSION 


Summary. We have studied in the present work sim¬ 
ple isotropic solids formed by permanent spring networks. 
Some relevant static properties have been reviewed and 
characterized in Secs. II A[ |IV| and |V A[ More impor¬ 
tantly, we have reconsidered theoretically (Sec. 0 and 
numerically (Sec. 0 the linear-response relation between 
the shear stress relaxation modulus G{t) and the shear 
stres s aut ocorrelation function G{t) both in the time 
(Sec. VB , Figs. ffl and the frequency domain (Sec. |V ^ 
Figs. lQ|andpT]). According to our key relations, Eqs. (0 


I 


G(tj and Cqt) must become different in the solid limit 
Geq > 0) and it is thus impossible to determine Geq 
from G{t) or its Fourier transforms G'{uj) and G"{uj). 
The short-time behavior of G{t) and the large-u; limit of 


G'{uj) only yield the stress-fluctuation contribution crp 
to the equilibrium modulus Geq = /^a — ctf (Sec. II A) as 
already stressed elsewhere mm- The peak Regen¬ 
cies observed for G' {uj) and G"{uj) (Figs, and 11) 
are qualitatively expected P^ITi] from the breakdown 
of continuum mechanics for large wavevectors q ^ 1 
seen from the Fourier transformed transverse displace¬ 
ment field (Fig. [^. 

Discussion. It is obviously common and may often 
be helpful to describe a plateau of G (t) at short or inter¬ 
mediate times (or equivalently at large or intermediate 
frequencies for G'(co’)) in terms of a finite shear modulus 
Gm of a dynamical relaxation model, such as the Maxwell 
model for viscoelastic fluids or the reptation model of 
entangled polymer melts mm- However, such a model 
allowing the theoretical interpretation of the data should 
not be confused with the proper measurement procedure^ 
Eq. 0. and a finite model parameter Gm not with the 
equilibrium modulus Geq of the system which, inciden¬ 
tally, vanishes both for a Maxwell fluid or a linear poly¬ 
mer melt. In this sense different operational “static” and 
“dynamical” definitions of the shear modulus are used in 
the literature for describing glass-forming liquids close 
to the glass transition [TTHini m- This may explain 
why qualitatively different temperature dependences — 
cusp singularity pn [22] vs. finite jump [fll [181 [SO] — 
have been predicted recently. Hence, while our recent at¬ 
tempts to determine Geq(T) for two glass-forming model 
systems [21] are consistent with a continuous cusp, this is 
not necessarily in contradiction with a jump singularity 
for Gm(T) determined from an intermediate shoulder of 

G{t) [23ES]. 

Outlook. We note finally that generalizing our key re¬ 
lations one obtains readily that 


M{t) = G{t) + Meq for t > 0 (48) 

for the relaxation modulus M{t) of any continuous inten¬ 
sive variable I with Meq = dl/dX being the equilibrium 
modulus and G{t) = f3V{6I{t)6I{0)) the corresponding 
autocorrelation function computed at a constraint ther¬ 
modynamically conjugated extensive variable X. In the 
frequency domain this leads to M'{uj) = Meq + G'{uj) 
and M"{uj) = G"{uj). A natural example is provided by 
the relaxation modulus M{t) = K{t) associated to the 
normal pressure I = P. While G{t) and G'{uj) must ob¬ 
viously vanish in the static limit for, respectively, large 
times and small frequencies, K{t) and K'{uj) approach a 
finite compression modulus Meq = i^eq foi* stable (non- 
critical) thermodynamic systems. 
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